module ComputeResidual
	
implicit none
	
	contains
	
	real(kind=8) FUNCTION computeLargestRelResid(folderName,fileNum)
	!在使用选手的作品对算例中的所有线性方程组进行完求解并输出结果后，可调用本函数计算所有线性方程组的解的最大残差
	use BasicMathSubroutines
	use ReadData, only: readMatrix, readVector
	
    IMPLICIT NONE
        character(len=100) :: folderName
		integer,intent(in) :: fileNum
		character(len=5) :: matrixSuffix = ".matr"
		character(len=4) :: rhsSuffix = ".rhs"
		character(len=10) :: solSuffix = ".submitSol"
		integer, parameter :: fileNamePrefixLength = 5
		character(len=100) :: formatDescription
		character(len=100) :: fileNamePrefix
		REAL(kind=8), DIMENSION(:,:), pointer :: matrix
		REAL(kind=8), DIMENSION(:), pointer :: sol,rhs
		integer :: i,j
		
		write(formatDescription,'(I1)') fileNamePrefixLength
		formatDescription = '(I'//TRIM(ADJUSTL(formatDescription))//')'
		computeLargestRelResid = 0.
		do i = 1, fileNum
			write(fileNamePrefix,formatDescription) i
			do j = 1, fileNamePrefixLength - LEN(TRIM(ADJUSTL(fileNamePrefix)))
				fileNamePrefix  = '0'//ADJUSTL(TRIM(fileNamePrefix))
			end do
			call readMatrix(TRIM(ADJUSTL(folderName)) // TRIM(ADJUSTL(fileNamePrefix)),matrixSuffix,matrix)
			call readVector(TRIM(ADJUSTL(folderName)) // TRIM(ADJUSTL(fileNamePrefix)),rhsSuffix,rhs)
			call readVector(TRIM(ADJUSTL(folderName)) // TRIM(ADJUSTL(fileNamePrefix)),solSuffix,sol)
			if (norm(matmul(matrix,sol) - rhs) / norm(rhs) > computeLargestRelResid) then
				computeLargestRelResid = norm(matmul(matrix,sol) - rhs) / norm(rhs)
			end if
		end do
		deallocate(matrix,sol,rhs)
		
        RETURN
	END FUNCTION computeLargestRelResid
	
	
	real(kind=8) FUNCTION computeLargestRelResidOfMultiplier(folderName,fileNum)
	!检查矩阵向量乘的精度
	use BasicMathSubroutines
	!use BLAS95
	use ReadData, only: readMatrix, readVector
	
    IMPLICIT NONE
        character(len=100) :: folderName
		integer,intent(in) :: fileNum
		character(len=5) :: matrixSuffix = ".matr"
		character(len=4) :: rhsSuffix = ".rhs"
		character(len=10) :: solSuffix = ".sol"
		integer, parameter :: fileNamePrefixLength = 5
		character(len=100) :: formatDescription
		character(len=100) :: fileNamePrefix
		REAL(kind=8), DIMENSION(:,:), pointer :: matrix
		REAL(kind=8), DIMENSION(:), pointer :: sol,rhsUserRoutine,rhsMatMul
		integer :: i,j
		integer :: matrixSize
		real(kind=8) :: tempResidual

		real(kind=8), dimension(:,:), allocatable :: localMatrix
		real(kind=8), dimension(:), allocatable :: localSol
		real(kind=8), dimension(:), allocatable :: localRhsUserRoutine
	
		write(formatDescription,'(I1)') fileNamePrefixLength
		formatDescription = '(I'//TRIM(ADJUSTL(formatDescription))//')'
		computeLargestRelResidOfMultiplier = 0.

		print*, "in computeLargestRelResidOfMultiplier"

		do i = 1, fileNum
			write(fileNamePrefix,formatDescription) i
			do j = 1, fileNamePrefixLength - LEN(TRIM(ADJUSTL(fileNamePrefix)))
				fileNamePrefix  = '0'//ADJUSTL(TRIM(fileNamePrefix))
			end do
			call readMatrix(TRIM(ADJUSTL(folderName)) // TRIM(ADJUSTL(fileNamePrefix)),matrixSuffix,matrix)
			call readVector(TRIM(ADJUSTL(folderName)) // TRIM(ADJUSTL(fileNamePrefix)),solSuffix,sol)
			
			allocate(rhsUserRoutine(size(sol,1)), rhsMatMul(size(sol,1)))
			
			matrixSize = size(sol,1)
			!print*, "mulMatrixVector round : ", i, ", fileNum : ", fileNum
			!print*, "sol size : ", size(sol,1), ", matrix size : ", size(matrix,1), ";",size(matrix,2), ", y size : ",  size(rhsUserRoutine, 1)
			! print*, "sol value : ", sol(1)

			allocate(localMatrix(matrixSize,matrixSize), localSol(matrixSize), localRhsUserRoutine(matrixSize))
			localMatrix(:,:) = matrix(:,:)
			localSol(:) = sol(:)
			localRhsUserRoutine(:) = rhsUserRoutine(:)

			!call gemv(matrix,sol,rhsUserRoutine)!选手需要替换此矩阵向量乘函数为作品
			call mulMatrixVector(matrixSize, localMatrix, matrixSize, localSol, 1, matrixSize, localRhsUserRoutine)!选手需要替换此矩阵向量乘函数为作品
			
			rhsUserRoutine(:) = localRhsUserRoutine(:)

			rhsMatMul = matmul(matrix,sol)
			tempResidual = norm(rhsMatMul - rhsUserRoutine) / norm(rhsMatMul)
			if (tempResidual > computeLargestRelResidOfMultiplier) then
				print*, "tempResidual : ", tempResidual
				computeLargestRelResidOfMultiplier = tempResidual
			end if

			deallocate (localMatrix,localSol,localRhsUserRoutine)
		end do
		
		deallocate(matrix,sol,rhsUserRoutine,rhsMatMul)
        RETURN
    END FUNCTION computeLargestRelResidOfMultiplier
	
	
end module